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Abstract 

We extend the POMWIG Monte Carlo generator developed by B. Cox and J. Forshaw, to 
include new models of central production through inclusive and exclusive Double Pomeron 
Exchange in proton-proton collisions. Double photon Exchange processes are described as 
well, both in proton-proton and heavy-ion collisions. In all contexts, various models have 
been implemented, allowing for comparisons and uncertainty evaluation and enabling detailed 
experimental simulations. 
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Program summary 

Title of the program: DPEMC, version 2 . 3 

Computer: any computer with the FORTRAN 77 compiler under the UNIX or Linux operating 
systems. 

Operating system: UNIX; Linux versions 7.x, 8.x 
Programming language used: FORTRAN 77 
High speed storage required: < 25 MB 

Keywords: Proton- (anti)proton collisions, diffraction, double pomeron exchange, double photon 
exchange. 

Nature of the physical problem: Proton diffraction at hadron colliders can manifest itself in many 
forms, and a variety of models exist that attempt to describe it 0121 113 IU 113 ■ This program im- 
plements some of the more significant ones, enabling the simulation of central particle production 
through color singlet exchange between interacting protons or antiprotons. 

Method of solution: The Monte-Carlo method is used to simulate all elementary 2^2 and 2 — >• 1 
processes available in HERWIG. The color singlet exchages implemented in DPEMC are implemented 
as functions reweighting the photon flux already present in HERWIG. 

Restriction on the complexity of the problem: The program relying extensively on HERWIG, the 
limitations are the same as in jH]. 

Typical running time: Approximate times on a 800 MHz Pentium III : 5-20 minutes per 10000 
unweighted events, depending on the process under consideration. 
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1 Introduction 



Doubly diffractive Higgs boson production has been much debated during the last decade, in a 
variety of theoretical approaches and predictions |TJ 12 E3 HI El • Now that the Tevatron is running, 
with both the DO and CDF experiments well equipped for forward proton detection, there is 
good hope to clarify the situation in the mid-term. Meanwhile, potential upgrades of the LHC 
experiments ATLAS and CMS are being discussed, and a systematic evaluation of the DPE physics 
potential becomes necessary The present work is meant as a step in this direction. 

Available programs for diffraction simulation are POMPYT |6j and RAPGAP jjj; the factorized 
model of Cox and Forshaw, POMWIG [4J; and the Soft Color Interaction model of Enberg, Ingelman 
and Timneanu [3]. POMWIG is implemented as an extension to HERWIG [S] and relies on the following 
branchings: 

p^p' + P, followed by P -> q + X (1) 

where p and p' are the incoming and outgoing protons, the Pomeron is denoted P, and q is the 
parton (quark or gluon) involved in the hard subprocess. In this inclusive picture, the final state 
can be made in lowest order of e.g. jets, photon pairs, lepton pairs, or a Higgs boson, accompanied 
by remnants from the IPP collision and the outgoing protons. 

Technically, the simulation of 77 interactions in ee collisions, provided by HERWIG, is trans- 
formed into pp induced PP interactions through a reweighting of the photon flux and a proper 
choice of the structure functions. The original implementation of this procedure, POMWIG, was 
based on a Regge fit to the Pomeron flux and Pomeron parton distribution functions provided by 
the HI experiment 0. 

Note that nothing prevents us from applying the same machinery to simulate exclusive PP 
interactions, i.e. processes where the entire PP center-of-mass energy is used in the production 
of a heavy central system. Most often, the Pomeron is seen as a gluonic object, so that the 
branchings 

p^p' + W, followed by P -> g + g NP (2) 

result in a hard gluon-gluon interaction, with e.g. jet pairs, photon pairs or Higgs bosons in the 
final state. The color screening soft gluon, g^p, couples to its counterpart from the opposite 
beam. Pomeron remnants are absent, and the protons again leave intact. Contrarily to the 
inclusive processes, which are being measured at the Tevatron (through dijet production |10j). 
the exclusive DPE has not yet been observed. This category of events is most appealing however, 
because Higgs boson production through this channel provides a very competitive measurement 
of its mass [TT], and possible determination of its spin and parity, which is a major challenge in 
non- diffractive channels. 

It is also simple to describe 77 interactions in the same framework : this is just going back to 
the original HERWIG where, to accommodate proton- or heavy ion-induced photon fluxes, we need 
to take into account the corresponding form factors. Compared to ee induced 77 interactions, 
an additional complication comes from the fact that one needs to simultaneously forbid strong 
interactions between the nuclei. To achieve this, one generally works in impact parameter space 
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a) b) c) 

Figure 1: a) Exclusive QCD, b) Inclusive QCD and c) QED mediated dim-active reaction in pp 
or Heavy-Ion collisions. 

and requires this impact parameter to be large enough for the nuclei not to geometrically overlap. 
These interactions have interest in their own right (they are a basic ingredient of the study of 
ultra-peripheral heavy- ion collisions, for example), but also provide a natural lower bound on the 
exclusive central production in pp collisions. In the worst case, these channels thus provide a 
conservative estimate of the physics potential of DPE. 

This flexible framework allows to implement a number of models alternative to that of Cox 
and Forshaw in an analogous way. We focus here on the exclusive model of Bialas and Landshoff 
(BL) P, and on its inclusive extension by Boonekamp, Peschanski and Royon (BPR) [2]. On 
the 77 side, we propose an implementation of the computations by Cahn and Jackson (CJ) [12], 
Drees, Ellis and Zeppenfeld (DEZ) [13], Papageorgiu [Tlj . and Budnev et al [THj . 

We can unfortunately not pretend to be exhaustive, and the implementation of further models 
is left for future work. One notable example is the model of Khoze, Martin and Ryskin j^J. We 
also escape discussing proton non-survival through soft parton rescattering; the cross-sections 
presented here should thus be understood as incomplete in that respect. Note however that the 
results of the inclusive DPE models can be rescaled to the measured DPE dijet cross-section, as in 
[2]. On the other hand, in QED mediated diffraction, the proton collisions occur at much higher 
impact parameter, and the probability of parton rescattering is small [3]. 

The next section gives an overview of the models newly implemented (we leave most physical 
discussions aside and orient the reader to the references quoted above). Some results are then 
given in the third section, and the Appendix provides a guide to the package, and lists all settings 
and switches that control the output. 

2 DPE models implemented in DPEMC 

This section describes the new models and processes which can be studied using DPEMC. The model 
of Cox and Forshaw is described in 0] and of course still available. 
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Central system 


Partonic process 


Higgs bosons 
Dijets 

Photon pairs 


gg^H 

gg -> gg, gg w 
gg -> 77 



Table 1: Hard subprocesses relevant to exclusive DPE. 



2.1 Exclusive DPE in the Bialas-Landshoff model [lj 

The production cross section for a heavy system X (such as a Higgs boson, a jet pair, or a photon 
pair) in exclusive DPE collisions: 

p + p — >p + X + p, (3) 
as illustrated in Fig. [T^i, takes the form [U |2] : 

At (cxc) = / «2 Fp/^i)Fp/ p (6) d*G70 - (4) 

The convolution in Eq. (jlj) is over the momentum fractions of the initial protons, £i and £2, 
carried away by the Pomerons. The flux of the Pomeron inside the proton is parametrized as: 



v 



2 



Fv/M^^Jdv 2 -^^, (5) 

where v 2 is related to the momentum transfer squared and the proton momentum loss through 
— t/(l — £). Here the parameters Af 2 = 288/7T 5 , (3 = 4GeV~ 2 and the Regge trajectory 



ajp(t) = «p(0) + apt, (6) 

with intercept ap(0) = 1.08 and slope ct F = 0.25 GeV -2 , are taken as in the original Bialas- 
Landshoff model, Ref. pQ. 

The partonic cross section da in Eq. (jlj) describes the possible hard subprocesses gg — > X, of 
which a list is given in Tabled in this context, the cross-sections satisfy color and angular mo- 
mentum selection rules; namely the initial state is color singlet, and only the helicity configuration 
with J z = and P = +1 contributes to a [HI]. These features have important consequences on the 
sensitivity to Higgs boson searches, since they imply important suppression of the backgrounds 
while being harmless to a naturally scalar and color-singlet signal. 

Note that properly speaking, the Bialas-Landshoff expressions have been derived for Higgs 
bosons and massive quarks only. We propose here a generalization to gluon pair production, by 
using the corresponding color- singlet, J z = subprocess cross-sections instead of the original ones. 



2.2 Inclusive DPE 

We implement the model of BPR as an alternative to that of Cox and Forshaw (CF) jl]. Although 
different in spirit, both lead to similar expressions and final states. The production cross-section 
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Central system 


Partonic process 


Higgs bosons 
Dijets 

Photon pairs 
Lepton pairs 


99, qq^ H 

99 -> gg(qq), qg qg, qq -> 00(5?), 99' -> 99' 

99 77, 99 -> 77 
g9 - ^ W 



Table 2: Hard subprocesses relevant to inclusive DPE. 

for a system X {X = Higgs, Dijet, Diphoton, Dilepton) in inclusive DPE collision (see Fig. Efc): 

p + p — > p + X + p + Pomeron remnants, (7) 
can be expressed as follows 0: 

da {inc) = ^2 dx i dx * d €i d & Fjp/pi^Fjp/pi^) fi/jp(x u /j 2 )f j/]P (x 2 ,fi 2 ) da(ij -> X). (8) 

i,3 

Besides the convolution over the proton momentum fractions, the integration runs over the 
Pomeron parton distribution functions. The Pomerons momentum fractions carried by the collid- 
ing partons are denoted X\ and x 2 - The proton induced Pomeron flux, is the same as in 
the exclusive case (the normalization and the parameters /3, ap(0) and ajp are assumed identical 
to those used in PQ). The partonic content of the Pomeron is expressed in terms of distribution 
functions: 

fi/ip(xi, fx 2 ) = Xi Gi/-p(xi, fx 2 ), (9) 

where Gi/-p(xi, fx 2 ) are the "true densities" of the partons inside the Pomeron as measured by the 
HI collaboration jU] . They are evaluated at a scale fx 2 given by the mass of the central system. The 
integral of fi/ip{xi, fx 2 ) is normalized to 1, so that in the limit fi/ip{xi, fx 2 ) — > S(xi) the exclusive 
expressions are recovered. 

Note that the cross-sections obtained in this way are rescaled in [2] according to the double 
diffractive dijet cross-section measured in ^Hj; this is not done in the program. 

Some of the possible hard subprocesses gg, qq — ► X are listed in Table El 

2.3 Two-photon reactions in pp and heavy- ion collisions [1121 113L flU I15| 

The production cross section for a central system X in heavy-ion collisions through double photon 
exchange (see Fig. Ht): 

N + N^N + X + N, (10) 
can be approximated by a factorized form: 
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Central system 


Partonic process 


Higgs 
Dijet 

Lepton pairs 
Photon pairs 


77 — > H 
77 -> qq 
77 — > 11 

77 -> 77 



Table 3: Hard subprocesses occurring in double photon exchange. 



J d&db F 7/N (6)F 7/N (6) ^(77 - X). 



The photon flux from an ion of electrical charge Z and mass Mjy is parametrized as follows by 
DEZ |I3| : 



aZ 2 1 



-* 2 M*/Q§ 



+ 1 



en 



Ei^/Q^ 



(12) 



where the momentum transfer scale Qq ~ 60MeV and Ei(x) = f°°dte t /t. Formula is 
valid for heavy nuclei only. An alternative formulation, derived in impact parameter space, yields 



2aZ 2 



2 Uo 



2, e 



6 



K 2 



6 



(13) 



In Eq. (Unj), i£o,i(£/£o) are the Bessel functions, while £o = 1/{RMn) with R the radius of the 
ion. As discussed in the Introduction, this parameter appears because the impact parameter space 
integration is limited to the region where the colliding nuclei do not geometrically overlap. One 
can use R ~ 1.2A 1 / 3 for collisions of ions of atomic number A, or the proton radius to estimate the 
77 cross-sections in proton-proton collisions. The partonic cross sections relevant in this context 
are listed in Table El 



3 Results and Perspectives 

All implemented fluxes are displayed in Figure El Apart from normalization, the CF and BL 
Pomeron fluxes differ mainly because of different choices of the Pomeron intercepts (ap(0) = 1.2 
is taken by CF, while «p(0) = 1.08 is used by BL and BPR). The nucleus induced photon flux is 
illustrated for the proton (f y / p ), Calcium (/ 7 /c a ), and Lead (/ 7 /p&). 

Figure El illustrates Higgs boson production cross-sections in various models, namely the in- 
clusive models of CF and BPR, the exclusive cross-section by BL, and the photon-mediated 
cross-section by Papageorgiu. Figure 0] displays the "background" in the BPR model (the total 
and bb cross-sections are displayed) and in the exclusive BL model. The exclusive bb prediction 
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Cox-Forshaw Pomeron 
Cox-Forshaw Reggeon 
Bialas-Landsholf Pomeron 




10 r 



10 =- 




0.0004 0.0006 



0.0008 0.001 




0.0008 0.001 




0.003 0.004 



Figure 2: Beam induced fluxes implemented in the program. Upper left : Pomeron flux according 
to [1] and PUIIIj and the Reggeon flux of [3]. Upper right : proton induced photon flux according 
to and [T3|. Lower left and right : Calcium and Lead induced photon flux following J2| and 

m- 

is much suppressed as expected, and is more distributed towards small angles than its inclusive 
counterpart. 

Considering the above examples and others that can be worked in an analogous way, it is 
possible using this program to perform complete studies of experimental sensitivity to diffractively 
produced Higgs bosons. These studies can be done following models of inclusive or exclusive 
double Pomeron exchange, and of double photon exchange in proton or heavy ion collisions. 
Supersymmetric enhancements of the quark and gluon couplings to the Higgs boson can also be 
included as in HERWIG. 

The implementation of further models of hard diffraction, such as that of Ref . j3] , will follow in 
subsequent versions of the program. Another important addition will be the inclusion of diffractive 
proton dissociation. Finally, since the beam particles currently need to be defined as electrons, it 
is impossible to simulate soft underlying events or multiple parton interactions. This shortcoming 
will also be handled in the near future. 

In summary, the aim of this work is not to promote one model or the other, but to provide a 
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CF (inc.) 
BPR (inc.) 
BL (exc.) 

Papageorgiu (exc.) 




50 100 150 200 250 300 350 

m H [GeV] 

Figure 3: Inclusive Higgs boson production cross-sections following CF |I] and BPR [2j. The 
exclusive predictions of BL 1] and Papageorgiu [Hj are also illustrated. 

tool that should allow to determine those that are compatible with the available and forthcoming 
data, and to prepare the analyses and constrain the predictions for the LHC. 
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A Manual 

The source code for the program can be obtained from the authors or from Ref. |17j . and in- 
cludes an example main program. Some aspects of the implementation of DPEMC are described 
in Appendix [Bj here we discuss how to run the program. The processes are controlled using the 
following flags: 

• IPR0C,ID,IL select the hard process as in HERWIG; 

• TYPEPR controls how the hard cross-section is to be evaluated; inclusively (TYPEPR^ INC) 
or only the color-singlet Jz = amplitude (TYPEPR= ' EXC ' ) ; 
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Figure 4: Distributions of the average dijet rapidity (yjj = {y\ + t/ 2 )/2) in the inclusive model of 
BPR compared to the exclusive model of BL [Tj. The total inclusive dijet rate is displayed, as 
well as its bb component, to be compared to the exclusive bb cross-section. 

• TYPINT switches between photon initiated processes (TYPINT='QED'), and gluon/quark ini- 
tiated processes (TYPINT= ' QCD ' ) ; 

• NSTRU controls the P or 7 flux; it can take values from 9 to 15, with meanings detailed in 
Table |3J For Pomeron fluxes, it controls the parton density functions to be used as well; 

• when NSTRU=11, GAPSPR and CDFFAC set the absolute cross-section normalizations. GAPSPR 
controls the gap survival probability for exclusive DPE processes (TYPEPR^EXC); typical 
values are 0.1 at the Tevatron, and 0.03 at the LHC [S]. CDFFAC controls the normalization 
of inclusive DPE processes (TYPEPR =;| INC), by comparison of the raw dijet computation to 
the measurement of ^U]; the value 3.8 was found in [2|. 

• when appropriate, AI0N and ZION set the atomic and proton number of the colliding nuclei. 

As in POMWIG, the beam particles need to be set as electrons in the main program, and the 
HERWIG facility generating 77 interactions is used and reweighted to the user's needs, depending 
on the above settings. 

Since NSTRU controls which flux is to be used in the cross-section evaluations, the introduction 
of TYPEPR and TYPINT may seem somewhat redundant. This is however needed for technical 
reasons. First, exclusive IP IP interactions can not be treated as such since the gluon initiating 
the hard process carries the full momentum of its parent "beam particle", provoking errors in 
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NSTRU 


Flux and structure functions 


9 


Factorized model, Pomeron flux 0] 


10 


Factorized model, Reggeon flux [3] 


11 


Bialas-Landshoff Pomeron flux QjJ |2| 


12 


QED flux from Calm, Jackson [TJ]; R ~ 1.2A 1 / 3 


13 


QED flux from Drees et al., valid for heavy ions only 13] 


14 


QED flux in pp collisions, from Papageorgiu ^3] 


15 


QED flux in pp collisions, from Budnev et al [HD 



Table 4: Possible values of NSTRU and their meaning. 

HERWIG. To circumvent this problem, we needed to modify a few HERWIG routines. For example, 
the routine evaluating 77 — > // has been modified so that it can compute the color-singlet, Jz = 
gg — > qq,gg cross-section, depending on the values of TYPINT and TYPEPR. 

For example, to simulate bb production in the exclusive Bialas-Landshoff model, one thus needs 
to request 77 — > bb via IPR0C=16005, set TYPINT= , QCD , ) TYPEPR= , EXC ) so that the color- singlet, 
Jz = gg — > qq cross-section is actually computed, and NSTRU=11 requesting the Bialas-Landshoff 
flux. An exhaustive list of consistent flag settings is given in Table |SJ These settings need to be 
followed strictly, the program is producing unpredictable results otherwise 1 . 

It was also needed to provide some additional routines, for the computation of the 77 — > H 
and color- singlet, Jz = gg — > 77 cross-sections. For completeness, a computation of 77 — > 77 
is also provided. A list of the new IPROC values, and of ID settings specific to DPEMC is given in 
Table [3 

Finally, we provide an additional routine that corrects the event record after each event is 
finalized (remember that all events are recorded as ee induced 77 interactions). This is needed for 
readability, but also to enable correct experimental simulations. When called, this routine resets 
the beam particles as protons, the secondary photon "beams" can be redefined as Pomerons or 
Reggeons, depending on the value of NSTRU, and the particles initiating the hard process, always 
stored as photons when TYPEPR= ' EXC ' , are redefined as gluons when TYPINT =, QCD'. 

B Structure of the program 

We give here a list of the FORTRAN files that come with DPEMC. They have been tested to comply 
with HERWIG version 6.5. The first set of files contain the original HERWIG routines, which needed 
modifications for our purpose: 

herwig6500 . f : the original source code of HERWIG; the routines HWSFUN, HWEGAM, HWEPRO and 
HWHQPM, superseded by the following files, are renamed to avoid conflicts; 

hwsf un65 . f : structure function evaluation; 

1 We have tried to make the program detect unsound settings as much as possible. 
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Final state 


IPROC 


TYPINT/TYPEPR 


ZI0N/AI0N 


NSTRU 






Inclusive DPE: 








Higgs bosons 


1600+ID 


QCD/INC 






9, 10, 11 


13 ij et s 


1500 


QCD/INC 






q in 11 


1 Dnf A Tl T~\ Q1TO 

LjcpLOIl pdllb 


1350+IL 


QCD/INC 






Q ID 11 

y, iu, ii 


U"ll fill r\Q lTC 


2200 


QCD/INC 






Q ID 11 

<J , IU, 11 






Exclusive DPE: 








Higgs bosons 


9900+ID 


QCD/EXC 






11 


Dijets 


6000+ID 


QCD/EXC 






11 


Photon pairs 


9800 


QCD/EXC 
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Double photon exchange: 






Higgs bosons 


9900+ID 


QED/EXC 


Z/A 


1/1 


12 or 13 ; 14 


Dijets 


6000+ID 


QED/EXC 


Z/A 


1/1 


12 or 13 ; 14 


Lepton pairs 


6006+IL 


QED/EXC 


Z/A 


1/1 


12 or 13 ; 14 


Photon pairs 


9800 


QED/EXC 


Z/A 


1/1 


12 or 13 ; 14 



Table 5: Processes available for simulation, and the corresponding settings. The values of IPROC, 
ID and IL are set as in HERWIG [Hj, except for the new settings detailed in the Table The 
flags should be set in the main program, of which examples can be obtained with the source code. 
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IPROC 


ID/IL 


Process 


1600+ID 


ID=1,...,6 


H — ► qq (resp. d,u,s,c,b,t) 




ID=7,8,9 


H — > l + l~ (resp. e + e~ ,fi + \i~ ,t + t~) 




ID=10,11 


H -> W + W~, ZZ 




ID=99 


All decay modes 


1350+IL 


IL=0 


All leptons (incl. neutrinos) 




IL=1,...,6 


resp. e,u e ,fi,u^,r,u T 



Table 6: Reminder of the ID.IL settings for Higgs boson and lepton pair production in HERWIG. 



IPROC 


ID/IL 


TYPINT 


Process 


6000+ID 


ID=1,...,6 


QED/QCD 


ll/gg -»• QQ (resp. d,u,s,c,b,t) 




ID=11 


QCD 


99 QQ (all flavours) 




ID=13 


QCD 


99 ~^ 99 + QQ ( & 11 flavours) 




ID=7,8,9 


QED 


77 — > Z + /~ (resp. e + e~,/z + /i~,T + T~) 


9900+ID 


as for IPR0C=1600 


QED/QCD 


77/#5 -»• # 


9800 




QED/QCD 


77M -»• 77 



Table 7: Special ID, TYPINT settings for IPR0C=6000, and the new processes IPR0C=9800 and 
9900; together with TYPEPR= ' EXC ' , these return the color-singlet Jz = cross-sections. 
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hwegam65 . f : chooses the flux and generates Pomeron or photon momentum fractions and virtu- 
alities; 

hwepro65 . f : decides on the hard process; 

hwhqpm65 . f : calculates the cross-section for gg — * qq, gg in the J z = mode, for 77 — » 
qq, W + W~ , as in the original implementation of HERWIG, and for 77 — » gg in the 
J 2 = mode. 

The following files contain the new routines needed to run DPEMC: 

flux65.f : implementation of the Pomeron and photon fluxes inside a proton or heavy ion; 

hlqcd.f : parton distribution functions inside the Pomeron; 

pomstr . f : user defined Pomeron structure function; 

hwhigp65 . f : calculates the cross sections for gg — > H and 77 — > H; 

hwhqpp . f : evaluates the cross sections for gg — > 77 in the case of exclusive diphoton production 
and for 77 — > 77 in the heavy-ion collisions; 

hwnorm . f : calculates the normalization of parton distribution functions inside the Pomeron; 

hwfxer.f : corrects the event record. 

In addition to the above mentioned files, a file vegas.f, containing the program VEGAS ^H] for 
numerical integration, is also needed to obtain the normalization of Pomeron parton distribution 
functions. The example main programs provided with the package use HBOOK for histogramming, 
and thus need to be linked against the CERN program library |19j . 

An example main program is in example. f, and produces 10000 exclusive cc events in the 
Bialas-Landshoff model. It is compiled by executing gmake -f makeExample. 
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